beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, tree = genome_tree)
genome_counts_filt <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts),]
rownames(genome_counts_filt) <- NULL
genome_gifts <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
dist <- genome_gifts %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
beta_q1f <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, dist = dist)
Permanova
#Richness
betadisper(beta_q0n$S, sample_metadata$region) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.09234 0.092344 14.312 999 0.001 ***
Residuals 56 0.36132 0.006452
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Daneborg Ittoqqortoormii
Daneborg 0.001
Ittoqqortoormii 0.00037875
adonis2(beta_q0n$S ~ region,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| region |
1 |
1.700798 |
0.2101443 |
14.89903 |
0.001 |
| Residual |
56 |
6.392678 |
0.7898557 |
NA |
NA |
| Total |
57 |
8.093476 |
1.0000000 |
NA |
NA |
#Neutral diversity
betadisper(beta_q1n$S, sample_metadata$region) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.54684 0.54684 36.556 999 0.001 ***
Residuals 56 0.83770 0.01496
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Daneborg Ittoqqortoormii
Daneborg 0.001
Ittoqqortoormii 1.2786e-07
adonis2(beta_q1n$S ~ region,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| region |
1 |
2.506002 |
0.261813 |
19.86154 |
0.001 |
| Residual |
56 |
7.065723 |
0.738187 |
NA |
NA |
| Total |
57 |
9.571725 |
1.000000 |
NA |
NA |
#Phylogenetic diversity
betadisper(beta_q1p$S, sample_metadata$region) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.15221 0.15221 13.724 999 0.001 ***
Residuals 56 0.62106 0.01109
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Daneborg Ittoqqortoormii
Daneborg 0.001
Ittoqqortoormii 0.00048622
adonis2(beta_q1p$S ~ region,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| region |
1 |
0.1736353 |
0.103117 |
6.438466 |
0.001 |
| Residual |
56 |
1.5102323 |
0.896883 |
NA |
NA |
| Total |
57 |
1.6838677 |
1.000000 |
NA |
NA |
#Functional diversity
betadisper(beta_q1f$S, sample_metadata$region) %>% permutest(., pairwise = TRUE)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000358 0.0003577 0.0791 999 0.786
Residuals 56 0.253218 0.0045218
Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
Daneborg Ittoqqortoormii
Daneborg 0.791
Ittoqqortoormii 0.77954
adonis2(beta_q1f$S ~ region,
data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$C))),
permutations = 999) %>%
broom::tidy() %>%
tt()
| term |
df |
SumOfSqs |
R2 |
statistic |
p.value |
| region |
1 |
0.03790714 |
0.08343941 |
5.097979 |
0.08 |
| Residual |
56 |
0.41640027 |
0.91656059 |
NA |
NA |
| Total |
57 |
0.45430741 |
1.00000000 |
NA |
NA |
Richness diversity plot
beta_q0n$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(region) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = region, fill = region)) +
scale_color_manual(name="Region",
values=c("#6A9AC3","#F3B942")) +
scale_fill_manual(name="Region",
values=c("#6A9AC350","#F3B94250")) +
geom_point(size = 4) +
# stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
)
Neutral diversity plot
beta_q1n$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(region) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = region, fill = region)) +
scale_color_manual(name="Region",
values=c("#6A9AC3","#F3B942")) +
scale_fill_manual(name="Region",
values=c("#6A9AC350","#F3B94250")) +
geom_point(size = 4) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
)
Phylogenetic diversity plot
beta_q1p$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(region) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = region, fill = region)) +
scale_color_manual(name="Region",
values=c("#6A9AC3","#F3B942")) +
scale_fill_manual(name="Region",
values=c("#6A9AC350","#F3B94250")) +
geom_point(size = 4) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
)
Functional diversity plot
beta_q1f$S %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(region) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = region, fill = region)) +
scale_color_manual(name="Region",
values=c("#6A9AC3","#F3B942")) +
scale_fill_manual(name="Region",
values=c("#6A9AC350","#F3B94250")) +
geom_point(size = 4) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9) +
theme_classic() +
theme(
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 12),
axis.title = element_text(size = 20, face = "bold"),
axis.text = element_text(face = "bold", size = 18),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 16),
legend.title = element_text(size = 18),
legend.position = "right", legend.box = "vertical"
)
